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FOREWORD 

This report summarizes the research initiated under NASA Research 
Grant No. NSG 1030 and continued under NASA Research Grant No. NSG 1111. 
The purpose of this work has been to develop global tropospheric pol- 
lution models that describe the transport and the physical and chemical 
processes occurring between the principal sources and sinks of CH^ and 
CO. 

The report is divided into two chapters describing (1) the model 
development and description, and (2) the results of long terra static 
chemical kinetic computer simulations and preliminary short term 
dynamic simulations. The work on the dynamic transport /chemistry 
model simulations is continuing. 

The author thanks the National Aeronautics and Space Administration 
for their financial support. In addition, the helpful discussions and 
assistance of Dr. Henry G. Reichle, Jr.’s group at the Langley Research 
Center are acknowledged. Specifically, Ms. Shirley Campbell provided 
invaluable assistance in the computer program development. Finally, 
acknowledgement is made to the National Center for Atmospheric Research, 
which is sponsored by the National Science Foundation, for computer 
time used in this research. 
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NOMENCLATURE 

3 

c molar concentration of species in the air, kmoles/meter 

3 ‘ 

molar concentration of species in the absorbing phase, kmoles/meter 

> 

eic, 

c* molar concentration of species in the air if the ocean and air phases 
were in equilibrium, kmoles/meter^ 

3 

C molar density of air, kmoles/meter 

2 

F molar convective flux, kmoles/meter /second 

2 

molar flux at boundary, kmoles/meter /second 

2 

J molar diffusive flux, kmoles/meter /second 

H Henry's Law constant, dimensionless 

k reaction rate constant - first order reaction, second 

second order reaction, meters^/kmole/second 

k|j^ mass transfer coefficient, meters/second 

-1 

K absorption coefficient, meter 

M Molecular weight, kgrams/kmole 

0( ) order of magnitude 
p pressure, newtons /meter^ 

r radial coordinate position, meters 

3 

R ideal gas law constant, 8.314 X 10 joules /kmole/°K 

3 

R molar species generation by chemical reaction, kmoles/meter /second 

2 

S species source at boundary, kmoles/meter /second 

t time, seconds 

T temperature, °K 

u velocity in the (j)-dlrection, meters /second 

velocity in the 6-direction, meters/second 


V 



V 


w velocity in the r-direction, meters /second 

X mole fraction (volumetric mixing ratio) of methane, dimensionless 

y mole fraction (volumetric mixing ratio) of formaldehyde, 

dimensionless 

z mole fraction (volumetric mixing ratio) of carbon monoxide, 

dlmens ionles s 

ai dimensionless constant. Equation (50) 

tt 2 dimensionless constant. Equation (51) 

A- incremental change in variable 

A representative grid interval, meters 

2 

e turbulent eddy diffusivity, meters /second 

2 3 

Go kinetic energy dissipation rate, meters /second 

0 coordinate position relative to equator, degrees 

4> coordinate position, degrees 

-1 

W vortlcity, second 

Subscripts 
G ground level 

1 grid point in 4>-coordinate 

j grid point in 6-coordinate 

k grid point in r-coordinate 

a absorbing phase 

P pollutant 

PS pollutant soil product 

S soil 

T tropopause 



Vi 


Superscripts 
n time level 

r r-direction 

V volumetric averaged 

X methane 

y formaldehyde 

z carbon monoxide 

0 9-direction 

(j)-dlrection 

- grid scale averaging operator 
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CHAPTER 1 

MODEL DEVELOPMENT AND DESCRIPTION 

Transport /chemistry models that describe the circulation of a species 
in the atmosphere can be extremely beneficial in understanding the physical 
and chemical processes occurring between the sources and sinks of a pol- 
lutant. In this chapter, such a model for the methane-carbon monoxide 
system is discussed. The model considers the physico-chemical action of 
these pollutants in the troposphere. This restriction is convenient 
since the tropopause provides a natural boundary across which little 
transport occurs. The data on sources and sinks for these pollutants 
is discussed, and the estimates of these strengths are based on the best 
available information relative to the major anthropogenic and natural 
contributions. The incorporation of the source-sink descriptions into 
the model is discussed in detail. 

The distribution and concentrations of methane and carbon monoxide 
in the atmosphere are interrelated by the chemical reactions in which 
they participate. A chemical kinetic model based on the pseudo-steady 
state approximation for the intermediate species and for inclusion in 
the species continuity equation was developed to account for these 
reactions. Therefore, mass conservation equations are only required 
for the methane and carbon monoxide. 

The numerical procedure employed to mathematically describe the 
transport/chemistry is a mass conservative scheme employing an integral 
flux approach. It is fourth-order accurate in space which is desirable 
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in simulating convective processes in three space dimensions. Since 
computer storage places restrictions on the scale of transport processes 
that are explicity calculated, smaller scale mixing is described using 
an artificial diffusivity. This is analogous to the concept of the 
artificial viscosity which is useful in the global circulation models. 

The results of the computer simulations using this model will be dis- 
cussed in the next chapter. 

INTRODUCTION 

Global transport/chemistry models of pollutants can be employed to 
analyze the circulation of pollution from its sources to sinks. Further- 
more, these models can be significant in placing anthropogenic sources 
in proper perspective on a global scale. In this chapter and the sub- 
sequent chapter, the development and the computer simulation results 
of a tropospheric global model for methane (CH^) and carbon monoxide 
(CO) are presented. 

The analysis is accomplished by geographically distributing the 
sources and sinks of CO, CH^, and CH^O, and simulating their convec- 
tive and diffusive transport by numerical solution on the computer 
of the three dimensional turbulent diffusion equation. The atmospheric 
phenomena of these species are coupled through atmospheric chemical 
reactions that occur. Thus, the species must be considered simulta- 
neously. Generally speaking, the oxidation of methane produces 
formaldehyde which decomposes to carbon monoxide. Other sources and 
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sinks of these pollutants are, of course, operating. 

Most analyses to the present have utilized a global residence time 
approach^^^* ’ Models incorporating a multiplicity of sources 

and sinks have not generally been attempted. The model that is described 
employs known source and sink strength data, the atmospheric chemistry 
of the pollutants in question, monthly averaged climatological data, 
and the turbulent diffusion equation for CH^ and CO to establish global 
concentration distributions. 

PHYSICO-CHEMICAL CONSIDEEATIONS 

The model development is restricted to'the troposphere. This is 
a logical boundary since the tropopause provides a natural surface thru 
which the rate of mass transfer is relatively low. Furthermore, the 
photolytlc decomposition of CO 2 appears to be unimportant as a source 
of CO in the troposphere and this enables one to decouple the CO 
transport from the CO 2 transport. The sources and sinks of methane, 
formaldehyde, and carbon monoxide will be briefly reviewed to provide 
a better appreciation of this complex system. Figure 1 illustrates 
the principal interactions that occur. 

Sources and Sinks of CH, 

4 

The sources and sinks of methane appear to be reasonably well 
understood at the present. The anthropogenic sources are largely the 
result of internal combustion engines and oil drilling and refinery 
operations. These emissions can be fairly well mapped based on auto- 

i 

mobile density and industrial activities. 



Stratosphere 



Figure 1 


Schematic of the Methane-Carbon Monoxide Interactions 
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The natural sources apparently far out strip the man made sources - 

the principal ones being decaying vegetation and other biological action. 

Some of this biological action occurs within marine environments, and 

as a result the surface waters of the oceans, bays, and rivers appear 

to be supersaturated with methane. Lamontagne, Swinnerton, Linnenbom, 

(7) 

and Smith have reported equivalent surface ocean and sea water 
concentrations about 1,2 to 1.7 times the corresponding atmospheric 
concentrations. Specifically in open tropical ocean waters, the surface 
concentrations (4.7 X 10 ^ ml/1) corresponded to an equilibrium atmos- 
pheric concentration of 1.80 ppm, whereas : the measured atmospheric con- 
centrations averaged 1.38 ppm. Bay and river waters appear to be even 
more heavily supersaturated. Their results specifically cited the 
following supersaturation ratios: Chesapeake Bay - 14.3, York River - 

21.2, Mississippi River - 5.67, Potomac River - 36.0. These values 
may also be affected by' local pollution problems. The data of Brooks 
and Sackett^ on the coastal waters of the Gulf of Mexico generally 
support Lamontagne et al’s results. However, they report that in the 
Yucutan area, where there is a major upwelling of deep water with low 
hydrocarbon concentration, the Gulf of Mexico acts as a sink for 
methane. 

The principal sink mechanism for methane appears to be in the 
homogeneous gas phase reaction of methane with hydroxyl radicals. 

CH^ + OH ^ CH^ + H 2 O 


( 1 ) 
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The methyl radical, can subsequently undergo reactions which result in 

formaldehyde .and ultimately in CO formation. Thus, this sink for. 

methane provides one natural source for formaldehyde and carbon monoxide, . 

The following sequence of reactions is responsible for producing the 

( 2 ) 

hydroxyl radical , 


0^ + h V -> 0(^D) + 0^ 

(2) 

0(^D) + H^O 20H 

(3) 

O(^) +M->-0 + M- 

(4) 

O+O^+M-J-O^ + M 

(5) 


It will be noted later that the hydroxyl and atomic oxygen (0) are also 
important in reactions with CO. 

Sources and Sinks of CH^P 

The anthropogenic sources of formaldehyde appear to be relatively 
small - the main ones being direct emission from automobile exhaust and 
formation during photochemical smog episodes. These estimates can be 
fairly reliably based on past auto' exhaust emission estimates and studies. 

The only apparent natural source for CH^O is from the methane oxida- 
tion just cited. Levy^^*^*^^ and McConnell, McElroy, and Wofsy^^^^ have 
suggested the following steps in the formation of formaldehyde by this 
mechanism. 


CH, + OH CH, + H^O 
4 3 2 


( 1 ) 
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CH^ + 0^ + M ^ CH^O^ + M 
CH2O2 + NO GH^O + NO2 
CH ^02 + CH2O2 2CH2O + O2 
CH^O + 02 '^ ^^ 2 ° ^ 


( 6 ) 

(7) 

(9) 


Thus, the formaldehyde formation and concentration is directly coupled 
to the methane distribution. 

The main sink for formaldehyde is in the photochemical decomposition 


and reaction with hydroxyl radicals. The following 


be importort 


(9,10,11) 


GH 2 O + h V CHO + H (10) 
CH 2 O + h V ^ H 2 + GO (11) 
GH 2 O + OH GHO + H 2 O (12) 


The production of GHO also leads to carbon monoxide formation via 


(13,14) 


CHO + O 2 CO + HO 2 


(13). 


Therefore, the source and sink distribution of formaldehyde is primarily 
due to homogeneous gas phase reactions and is coupled to the methane and 
carbon monoxide distributions. 
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Sources and Sinks of CO 

( 2 ) ( 12 ) 

Bortner, Kuimoler, and Jaffe and, more recently, Seiler have 

summarized the sources, sinks, and concentrations of carbon monoxide. 

A very significant feature of the global distribution of CO is the dif- 
ference of the mixing ratios found in the Northern and Southern Hemi- 
spheres. North of the intertropical convergence zone typical concen- 
trations over oceans are 0.15 ppm - 0.20 ppm, whereas south of the 

intertropical convergence zone the CO mixing ratios drop rather rapidly 

( 12 ) 

to 0.08 ppm. Seiler' ' has summarized most of the data on which these 
results are based, and he suggests average tropospheric concentrations 
of 0.15 ppm and 0,08 ppm in the Northern and Southern Hemispheres, 
respectivly. 

The world-wide anthropogenic sources are estimated to be in excess 
of 300 million tons/year with two-thirds or more resulting from motor 
vehicle emissions. The remainder is distributed between stationary 
combustion sources, industrial processing, and incineration. Therefore, 
these sources are distributed largely according to motor vehicle density. 

The major natural sources of carbon monoxide appear to be the 
oceans, forest fires, terpene photochemistry, and gas phase reactions. 
The studies of Junge, Seiler, and Warneck^^\ Seiler and Junge^^^\ 
Swinnerton, Linnenbom, and Gheck^^^^ , Lamontagne, Swinnerton, and 
Linnenbom^^^^ and Swinnerton and Lamontagne ' indicate that the level 
of excess CO in the ocean corresponds to an equilibrium air phase con- 
centration of about 3,5 ppm. These data were obtained during ocean 
cruises. The recent data of Meadows and Spedding^^^^ 


however , sugges t 
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that the degree of supersaturation may not be nearly as large. The 
previously mentioned studies based their calculations on CO solubility 
obtained for pure CO in the gas phase at pressures greater than or equal 
to one atmosphere. Meadows and Spedding, on the other hand, conducted 
experiments in distilled water and seawater where the gas phase CO con- 
centration was in the 3 ppm - 18 ppm range and found the solubilities to 
be more than eight times as high. These data would suggest that the 
oceanic source strength is not nearly as large as the other workers 
have suggested and that the equilibrium air phase concentration may 
be more like 0.3 ppm to 0.4 ppm rather than 3.5 ppm. 

' It is, of course, reasonable to expect that the river, lake, or 
ocean regions near urban areas, where the atmospheric CO concentrations 
may be considerably larger than 3.5 ppm, may act as a sink for CO. 
Furthermore, it is possible that the oceans at the high latitudes serve 
as a sink for CO produced by oceans at the low latitudes. This would be 
expected since the warmer tropical waters would likely have a higher 
biological activity producing more CO. With the warmer water, the 
solubility 'of CO is reduced. Transport over colder waters with greater 
CO solubility capacity would create the possibility of these sections 
acting as a sink for CO. Thus, it is plausible that the oceans act 
as both a source and sink for CO, Similar arguments could be made 
for methane. 

In addition to terpene photochemistry and forest fires (combined 
sources are estimated at 23 x 10^ ton/year^^^^), another principal 
natural source of CO appears to be the gas phase reactions cited 
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earlier. The estimate of the gas phase oxidation of CH^ to CO is 

directly related to the hydroxyl concentration. Since there have been 

no direct measurements made of the hydroxyl radical concentration, 

there is considerable uncertainty as to the magnitude of this source 

of CO, Various estimates place this anywhere up to ten times the 

(12 19) 

anthropogenic source * . However, the relatively slow inter- 

hemispheric transport from the north to south and the substantial 
differences in the CO concentration between the two hemispheres would 
suggest much lower estimates since the gas phase oxidation would not 
be favored in either hemisphere. The reactions forming CO cannot be 
divorced from the reactions which consimie CO, of which the following 
seem to be important 


CO + OH CO^ + H (14) 
CO + 0 + M CO 2 + M (15) 
CO + N^O CO 2 + (16) 
CO + HO^ ^ CO 2 + OH (17) 


Reaction is reportedly first order in CO but zeroth 

order in N 2 O. This is a surface catalyzed reaction and would require, 
for complete accuracy, detailed information on the atmospheric aerosol 
as to size distribution and chemical composition. Since these data 
are very scarce, there. is considerable uncertainty in employing 
Reaction (16) in any model. 
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(26) 

Reaction (17) has been suggested by Westenberg' ^ to be important 

in atmospheric pollution problems. However, the results of Davis, Wong, 
(27) 

Payne, and Stief indicate that it is unimportant in the overall 
oxidation processes of CO. As a result, the present model development 
ignores Reaction (17) . 

Another natural sink of CO of seemingly large significance is the 

soil. The work of Inman^^^\ Liebl^^^\ Inman, Ingersoll, and Levy^^*^\ 

(31) (37) 

Ingersoll and Inman , and Ingersoll, Inman , and Fisher “ , point 

up this significance. The field studies of Inman, Ingersoll and co- 

workers exposed soils in situ and in the laboratory to test atmospheres 

containing 100 ppm of CO. These showed average uptake rates that varied 

from 1.1 mg CO/hr m to 64.5 mg CO/hr m . On the low end were soils 

under cultivation and desert areas, and on the high end were tropical 

deciduous forest areas and soils near roadways. By using an average 

concentration driving force of 50 ppm CO and assuming conditions of 

atmospheric pressure and 20®C, the mass transfer coefficient at the 

“6 “4 

surface corresponds to 5.27 x 10 ra/sec and 3,09 x 10 m/sec, 

(29) 

respectively. The data of Liebl were obtained at 

more realistic CO mixing ratios - on the order of 0.25 ppm - than that 

of Inman, Ingersoll and co-workers. Liebl also showed uptake of CO by 

the soils. In addition, he showed a CO production capability when the 

air above the soil was initially free of any CO. These data indicate, 

(13) 

as Seiler and Junge have suggested, that at low concentrations 
(around 0.2 ppm at 25'’'C) a temperature dependent equilibrium of CO above 
soils occurs. However, for soil temperatures below 20‘‘C, the CO mixing 
ratio continued to decrease to a value lower than 2 ppb. If this is 
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truly the case, then the soils can act as either sinks or sources for 
carbon monoxide in much the same manner as do the oceans. Evaluation 
of exchange coefficients based on Liebl's data indicates a value an 
order of magnitude or so larger. However, this may simply be indicative 
of soil difference since he used soil obtained from a greenhouse. 

These scattered values are probably indicative of the rate of the 
biological reaction that is occurring near the surface of the soil and 
may simply be representative of the type and concentration of the soil 
micro-organisms utilizing CO. One of the big uncertainities relative 
to soil scavenging is the determination of which fraction of these 

( 33 ) 

micro-organisms are anaerobic methane-producing and what fraction 
are aerobic CO^ producing. 

A source and/or sink common to all three species is leakage from 

and/or to the stratosphere. Since the Initial model is restricted to 

the troposphere, this leakage must be considered as sources and/or 

sinks. In the specific case of CO, leakage from the troposphere to 

the stratosphere is reasonable to expect since the CO that escapes 

thru the tropopause will typcially undergo chemical reactions and 

not return to the troposphere as CO. This is substantiated by the 

vertical profiles of CO which show a decrease of CO mixing ratio 

( 34 ) 

with height above the tropopause . The leakage of CH^ thru the 
tropopause appears to be similar to that for CO since vertical profiles 
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I 

MATHEMATICAL MODEL 

The transport of a gaseous compound in the atmosphere is mathe- 
matically described by the species continuity equation with the inclusion 
of any homogeneous generation or loss terms. By distributing the sources 
and sinks of the various species on the Earth's surface and at the tropo- 
pause consistent with the physico-chemical considerations, the appropriate 
boundary conditions can be Incorporated into the solution. 

Mathematical Overview 

The diffusion equation written in spherical coordinates using the 
turbulent eddy diffusivity concept is 


9x. 9x. V 9x. u 9x, 

A- , ^ ^ _i_ ^ 

9t 9r r 90 r cos 0 90 


E ,1 9 f 2„ ^^A, , 1 9 , „ „ 

C 2 9r ^ 9r ^ 2 ^ dQ ® ^ 96 ^ 

r r cos 0 




(18). 


In Equation (18) the eddy diffusivity has been assumed constant, and x^ 
is the mole fraction of species A, expressible as ppm by volume if so 
desired; C is the molar density, which for air can be determined by 
using the ideal gas law (C = -^) ; is the generation or loss of species 

/s 

A by chemical reaction; and is a term which accounts for the effects 

(37) 

of turbulence on the chemical reaction . The latter term arises 
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from turbulent fluctuations in the species concentrations. For a single 
first order chemical reaction R/ is zero, but for more complex chemical 
kinetics this term is not zero. However, in many cases this term is 
relatively unimportant and in the future development^ will be con- 
sidered to be zero. 

Due to the well known closure problem caused by the introduction of 
fluctuating quantities, one is usually forced to parameterize the turbulent 
quantities. Perhaps the easiest and most commonly used procedure to solve 
the closure problem is to introduce the eddy diffusivity as has been done 
in Equation (18) . There is one problem that should be recognized with 
the use of the eddy diffusivity. Since the time scale associated with 
atmospheric motion can be quite large, unlike that in the more frequently 
encountered boundary layer flows, the duration of time averaging must 
generally be relatively long to include all fluctuations. In fact, this 
length of time is prohibitively long, and one must be satisfied with 
averaging times that are short in comparison. As a result, the eddy 
diffusivities that are used must be consistent with this concept and 
are directly dependent upon the time and length scale of the averaging. 

As Equation (18) now stands, numerical solution is required. 

However, it is further complicated by the highly non-linear character 
of the generation term. The problems associated with the numer|.cal 
solution are discussed in the next section. 
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Chemical Reaction Model 

Based on the considerations presented previously, the following 
reaction sequence currently seems plausible to describe the generation 
terms. 


0^ + hv 0( d) + 0. 


0( d) + HgO 2 OH 


1 ^ 

0( d) + ^ 0 + 


0 + Og + 0^ + 


GH^ + OH CH^ + HgO 


GH^ + Og + Mg, ^ CH^ 0^+ Mg 


CH^ Og + NO ^ CH^ 0 + NOg 


2GH^ ©2 2CH^ 0 


GH^ 0 + Og ^ CH 2 O + HOg 


^10 

CH_0 + hv + GO + 
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^11 

CH^O + h V ^ CHO + H 


'"12 

CH 2 O + OH CHO + HgO 


CHO + Og ^ CO + HOg 


^14 

CO + 0 + ^ CO 2 + 


^15 

CO + OH ^ CO 2 + H 


he 


CO + N„0 ^ CO + N, 

^Surface ' 


By considering this mechanism, the generation term for methane can be 
written as 


^CHj^ °0H (19) 

and for CH^O it becomes 

^CH20 " ^9°CKjO°02 ~ ^^10 ^11^°CH20 

■ "12 °CH^O « 0 ) 

Finally, the generation term for CO is 

®G0 “ ^10 °CH20 '*■ ^13 '^CHO " ^14 ®CO °0 
” ^15 °C0 °OH ■ ’^16 °C0 


( 21 ) 
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In order to evaluate Equations (19) thru (21) , it is necessary to 
know the concentrations of the various radical species involved. The 
pseudo-steady state approximation can be used to obtain these estimates. 

Let us digress momentarily and discuss the pseudo-steady state approx- 
imation since it has only been Infrequently used for geochemical pro- 
blems. If one considers the reaction sequence to .be occurring in a 
static system, the set of ordinary differential equations that de- 
scribes the chemical kinetics is representative of a stiff system. 

A system of differential equations is said to be stiff if the charac- 
teristic time constants of the individual steps vary greatly, i.e. , 
the rate constants are substantially different. Such is the case with 
the current system. Characteristic of the solution behavior of these 
systems is a very rapid initial change in some variables (these are 
said to be the stiff variables) followed by a slowly varying state. 

Over this time period, the non-stiff variables may change little or 
not at all. 

The numerical solution of stiff systems presents great computational 
difficulties, and a substantial literature (c.f.. Reference (38) as a 
recent summation of the state of the art) has been developed around 
their solution. For absolute numerical stability of the system, most 
integration schemes require that the time step be less than the smallest 
characteristic time constant. Numerical integration schemes such as 

(39) 

the Gear package are attempts at circumventing this problem. 

The pseudo-steady state approximation is another way of circumventing 
this problem if the variables of interest are the non-stiff ones. For the 
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current situation, this is true since our principal interests are with 
CH^ and CO. By exploiting the pseudo-steady state approximation, we can 
reduce a stiff system to a non-stiff system with little change in the 
accuracy of the non-stiff variables. As an -aside^ Lapidus, Aiken, and 
Liu^^^^ have shown that the pseudo-steady state approximation can 
actually improve the accuracy in certain very stiff systems where the 
computation time, by practical considerations, limits the smallness of 
the integration time step. 

The intermediate species, or stiff variables, exist in very low 

concentrations due to thelf high reactivity. As such they adjust to 

perturbations in the concentrations of CH^, CH 2 O, or CO very rapidly. 

These transients only last for fractions of a second (on the order of 
“5 

10 seconds and smaller) and the pseudo-steady state approximation 
is thus valid for time scales larger than this. This will be discussed 
in more detail later. 

Using the pseudo-steady state approximation and assuming that the 
concentrations of H^O, 0^, and the third bodies (M^, M^, can be 

reliably specified, the results are 

(kio + k^j^) 


^12''cH20^ 

^S^CH^ ^12^CH 0 ^IS^CO^ ^^2^H20 


2 3 4 


‘5‘=CH, ’‘12=CH 0 *‘15‘=C0> > 


( 23 ), 
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and 


\o " ^^10 ^ ^11^''cH 20 ■'■ a^c^„ '+ 'k\c ~ ^ 


5"CH^ ' "12"-CH20 ■*■ ^15^^00^ 


J J 14 


-k,,C,„ (24). 


<‘'2‘=H,0 + '^3 % > <Vo, X ■" *'14 °C0 “ 

2 3 2 4 14 


Interaction with Oceans 

In solving the turbulent diffusion equation, the boundary condition 
establishes whether the surface acts as a source or sink, or is passive 
to the transport process. This boundary condition is written as 


- e 


3c t 


r 3r' 


r=R 




- c. 


(25). 


r=R 


In Equation (25), c is the air phase concentration of the species, c. 

\ 

is the bulk concentration of the species in the ocean phase, c„ | is 

^ r=R 

the concentration on the ocean phase side of the interface, and is the 
ocean phase turbulent exchange coefficient. For sufficiently dilute 
systems, Henry’s Law can be used to relate gas phase and liquid phase 
concentrations at the interface. If H is Henry’s Law constant, then 
for equilibrium at the ocean-air interface 


‘^lr=R “ ® ^f,lr=R 


(26) 
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Thus the boundary condition is written as 


9c 

9r 


r=R 








c„) 


(27), 


where Cj represents the air phase concentration if the ocean and air 

Aj 

were in equilibrium; i,e.. 




(28), 


From this point the absorption coefficient, K^, is used for convenience. 




(29) 


I * 

Absorption (ocean operating as a sink) would be determined by c ' r=R^^a’ 
and desorption (ocean operating as a source) would be determined by 

The estimation of presents some difficulties and uncertainties. 

Henry's Law constant may be the least uncertain quantity and is known 

fairly accurately for distribution between air and fresh water . H 

depends on the water temperature and generally increases for increasing 

water temperature. Obviously this effect must be taken into account. 

f4** ^ 

Less certain is the estimate of kjj^. Okubo presents representative 

values of the vertical oceanic eddy diffusivity in the Cape Kennedy 
2 2 

area to be 1.3 cm /sec to 10 cm /sec. If the film model for mass 
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trasnfer is used to estimate and if the film thickness is assumed to 

be the 70 meter thick upper-ocean layer, then using the mid-range that 

2 

Okubo has cited (i.e. , 10 cm /sec) 
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2 , 

cm / sec 


70 m 


= 1.43 X 10 ^ m/sec 


(30), 


This value of agrees reasonably well with values estimated from 

Broecker, Li, and Peng’s report of the CO 2 flux (k,^= 3.31 x 10 ^ 

m/sec) and Junge, Seiler, and Warneck’s^^^ estimate of the oceanic CO 

-5 

source strength (k^= 0.410 x 10 m/sec). Values of which must be 

based on the grid mesh size in order to realistically simulate the 

subgrid scale mixing processes will be discussed later. 

* 


As long as H and c^ are known for the particular species, the 

boundary condition then accounts for the operation of the ocean as 

a source or a sink. For example, for CO appears to be around 3.5 

(17) 

ppm (although the results of Meadows and Spedding indicate a 
possibly much lower value) and for CH^ it is around 1.8 ppm. 


Interaction with Soils 

Interaction with soils is a combination of physical, chemical, 
and biological actions. However, if one assumes a very simple scheme 
of a general pollutant P interacting with the soil S, a result analogous 
to that for the oceans is obtained. For example, consider 

k’ 

s 

P + S PS (31) 

k" 
s 

PS P + S 


(32). 
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Then 


9c 


-e 


3r' 


■ H' =PS - K =S=P 


r=R 


(33) . 


If Cg and Cpg are relatively constant and in large excess one obtains 


9c^ 


-e 


9r 


= (Cg - Cy) 


(34), 


r=R 


:k 

where c^ is the equilibrium concentration. This is the quantity to 

(13) 

which Seiler and Junge ascribe the value of about 0.2 ppm at 25°C 

for CO. Furthermore, corresponds to the mass transfer coefficient 

(31) (29) 

determined from the work of Ingersoll and Inman' and Liebl' with 
CO, Unfortunately, data determined with CO cannot be extended to other 
species, such as CH^, as is possible with the ocean absorption. This 
is readily appreciated since the soil’s chemical and biological action 
would be different for each species. 

Leakage to the Troposphere 

Leakage to or from the troposphere creates what may be classified 
as artificial sources or sinks. If a complete atmospheric global dis- 
persion model were the object, these effects would not be realized 
as sources and sinks. However, model complexity precludes such con- 
siderations in the initial simulation. 

As a first approximation, the tropopause can be considered as a 
zero-flux boundary. This would be expressed as 
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-e 


8c 

3r 


= F„ 


r=Tropopause 


(35), 


where is zero for the zero flux boundary condition. However, most 

investigations indicate that F^, although small, is not zero. Specifically, 

the analyses of Seiler and Warneck^^^^ and Machta^^^^ indicate that for CO, 
—11 2 

~ 5 X 10 gCO/m sec. With this value of F^ the concentration gradient 

is quite small and changes in the CO concentration in the upper kilometer 
or so of the troposphere are less than 0.1%, even for modest values of the 
eddy diffusivity. Thus, a zero flux condition is a valid first approx- 
imation. 


NUMERICAL MODEL 


In the previous section, the physico-chemical model was described-. 
In this section, the numerical treatment is discussed. 


Convective Difference Schemes 

One of the primary objectives of numerical schemes for convective 
problems is to preserve conservation of mass. Some procedures, there- 
fore, employ a flux approach to numerically simulate the problem rather 
than directly finite-differencing the partial differential equation. ’ ’ 

These methods have been found very useful in numerical solution of j 

the Navier-Stokes equations and have received significant attention!^^^ : 

f 

especially with regard to turbulent flow calculations and atmospheric- 

) I 

fluid motion problems. One of the advantages to solving Equation (1.8) , 
as compared to the Navier-Stokes equations, is that it has linear 


\ 

convective terms whereas 'the latter contains nonlinear convective ' 
terms. However, the generation terms in Equation (18) do introduce 
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non-linearities in the proposed studies. 

Consider the difference element illustrated in Figure 2. By 
applying the conservation of mass for any particular species the fol- 
lowing equation can be written. 


n+1 n 

c. - c, . , r, 

. x,j,k i,3.k. 2 

'• .At ^ ® k 


cos 0 , Ad) A0 Ar = (F 
1 


<j) 


+ J 


<j> 


A Aft 

-F-,t . J..T . ,) 4r, A6 Ar + (cos 0, - F. , ,+ cos 6. -J. . - , 

1+1,3, k k 3-1 i,3-l,k 3-I i,3-l,k 




'tn*i,j,k+l 'k+1 '^i,i,k+l^^ \,j,k ^’^k 0jM4SAr (36) 


In Equation (36), i, j, and k refer to the grid point in the (f), 0, and r 
directions, respectively, and n refers to the time level. It should be 
noted that Equation (36) has been applied over two grid intervals in each 
direction. f 1* , , , F® . , , and f 5 . , are the convective fluxes in the 
superscripted directions and can be related to the velocities in these 
directions. Following Roberts and Weiss these can be expressed as 


i,3,k 


1 

4r, A0ArAt 
k 


^n+1 ®3+l ’^k+l 

SIS ‘^(^i»0>’^>t)u(<!)^,e,r)r 


dr d0 dt 


(37), 
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n+1 ‘^’i+1 ^+1 


i,j,k ^ 4r“cosi7A?ArAt ^ / / / c((f>,6^ ,r,t)v(<i),0 . ,r)cos9r dr d(j. dt 

^ . t " '^i-1 Vl (38), 


and 


i»j»k 


n+1 ^i+1 ® 3 +l 


/ / 


4r, COS0 .A(j)A6At n 4>. i 9- i 


/ c((j),0 jrj^,t)w(4>j0,r^)r cos0 d0 d({) dt 


(39). 


Development of Equations (37) thru (39) in difference form is the central 

point in establishing the numerical accuracy level, and for the current 

work it is fourth order in the convective’ terms. , , , jt . , , and 

i,3,k’ 

0 

j represent the diffusive fluxes in each of the superscripted 
directions. These are expressible in gradient form. 

Anticipating the development of the convective fluxes, the diffusive 
terms will be evaluated using the intermediate time level, n + 1/2, as 
well as the time levels n and n+1. After differencing the diffusive flux 
terms. Equation (36) can be rewritten as 


V C, , ^ , ,, C, , ,, + (F*^ 


- F*^ ) 

“*^4 -LI 


A t 


ijj»k i,j,k i,j,k i,j,k ^ i-1, j ,k i+ljjjk'^ 2r, cos 0 . A4> 

K J 


+ (cos 0 . - F® . T - cos 9 - F® . - . ) - 7 : — Q 

j-1 i,j-l,k 3+1 i,3+l,k'^ 2r cos 6.A0 


+ (r^ - r^ F^ ^ A t 

1 ■*'4 4 -1 -^4 4 ^ 


k-1 i,j,k-l "k+1 x,j,k+l" 23-2 ^ ^ 

k 


, , n+1/2 . n+1/2 n+1 n ^^r^i.i.k^ 

+ (x, . , + X. X. . - - X. . - ) 

i>3»k~l x,j,k x ,3 ,k Ar^ 



27 


^ ^^i,j,k+l ^i,j,k-l^ r^ Ar 


+ (C C '> s ^ 

^'^i,j,kfi "^i,j,k-r ^ x,j,k+i ■ ''i,j,k-l^ ^2 


+ - x^^ - x^ ') "9"'i,j,k' 


6^C_. , At 


r, A0 
k 


^ (C. ,. - C. ,.) ,. - xfi/? „) "a"' 


i.j+l,k x,j-l,k^ ^ i,j+l,k i,j-l,k^ 2 

k 


_ , n+1/2 _ ^n+1/2 . ^9^i,j V^n+1/2 . n+1/2 

^^i,j+l,k ^ij~l,k^ ,^2 ■*■,• ^Vl,j,k "■ Vl,j, 


2r, A0 
k 


_i_i ^aC . . , At 

_ n+1 n . <|) x, 3 ,k . 

i,j,k .2_2^ ,^2 ^‘^i+l,j,k 


c, cos 0 , Ad) 
k 2 


Vl,j,k^ '' i+l,j,k 


n+1/2 , 

X. . , , ) 

1-1. J.k 


e. At 

Z±. 


4r,^cos^0 ,Ad) 
k 2 


y- + rY . , At 
2 1,1, k 


(40), 


In Equation (40) , xY . , represents the mole fraction of the species 

1,3 

being considered, and C. . , represents the molar density of the air 

1,3 ,x 

determined from the prevailing temperature and pressure. 

The formulation of the diffusive flux difference terms uses a 
central difference scheme, which is essentially second order accurate. 
Roberts and Weiss state that this method is unconditionally stable 
and sufficiently accurate as long as e is small. Certainly in the 
horizontal motions this is true; however, subgrid scale motions usually 
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will dominate the vertical fluxes, and may not be considered small. 
Thus, it could be qualitatively argued that the procedure is essentially 
fourth order accurate in 6 and (}> but only second order accurate in r, 
as well as time. However, this drawback is easily overcome since the 
mesh size in the radial or vertical direction will be considerably smaller 
than the mesh length in the horizontal. Therefore, it appears that the 
accuracy of the scheme is still limited by the order of accuracy of the 
horizontal motions . 

The generation term in Equation (40) is a volume average over the 
region i+1, j+l> h+l> and over time At. Thus, it can be expressed as 

^n+l*^i+l®j+l ’^hd-l 

rY . , = —9 f f I f R((j),9,r, t)r^cos0drd0d(j)dt (41). 

8 r COS0 A<j)A9ArAt ^ a ^ 

^ ^ t “ ^-1 

It should be noted at this point that the generation terms will be 
determined from concentrations already calculated. An alternative 
scheme would be to calculate the concentration of one species at the 
new time level using the intermediate time level. Then this updated con- 
centration could be employed in the generation term for the next species. 
This might improve the computational stability characteristics but would 
have to be established by numerical experiments. These types of schemes 
are necessary in order to decouple the three diffusion equations for each 
of the three species. Otherwise, the three equations would have to be 
solved simultaneously, which, due to their nature, would either require 
linearization or some Iterative scheme. It is felt that these procedures 
are not justified. The individual generation terms can then be found 
from Equation (22) thru Equation (24) . 
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It is now appropriate to further consider the flux and generation 
terms, thus far only expressed in integral form. By expanding the 
product of concentration and velocity in a Taylor series in space and 
time, any desired accuracy can be achieved. For each of the flux terms, 
the following results can he obtained. 


f1* 




= (cu)f 




+ 1/3 (|§ 1^) 


n+1/2 


39 30'^l,j,k 


A0^ + 1/3 (!§”), 


n+1/2 


^3r 3r i,j,k 




+ 0(At^, A0^"P) 


(42) 


F® = (cv)® +1/3 A({)^ 

i>l»k x,j,k 3(j) 3^ x,j,k 


+ 1/3 

'^3r 3r'’ 


n+1/2 


Ar^ + 


0(At^, 


A<|)^r^ 


(43) 


x,j,k 


F 


r 

i>3,h 


(cw) 


r A9 

i, j ,k sin A0 


„ ~ n+1/2 

+ 1/3 IH.) 

^ 3([)^i,j,k 


sin A6 
A0 



, „ ,3c 2-A0^ sin A0 ^ 

^ 2 % 30^x,j,k at— ^ 


+ 0(At^, A0P (j)^“P) 


(44) 


We will be content with fourth-roder accuracy in the space variable and 

2 

second order in time as indicated by terms 0(At , A-) . It should be 
noted that in order to' assure second order accuracy in time, the direction 
of integration must reverse with each time step. Furthermore, for 9.> 45' 
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the accuracy in the 0— space variable is decreased. However, at 0^ as 
large as 83° the accuracy is apparently still third-order. 

Roberts and Weiss suggest expressing the zeroth-order terms 

in a manner as illustrated below. 


(cu)^ . t = 1/6 [u. . . - c? . )] 


(45) 


Equation (45) alternates with the following form for the reverse direction 
integration. 


(cu)f . , = 1/6 [u. . (8c?^y^ - c? - c?;|;J- . . )] 


(46) 


The derivatives in Equations (42), (43), and (44) are differenced using 
a central difference scheme. 

The generation term described by Equation (41) can be handled in a 
manner slightly different than the flux terms. The result is 


TT rv 4 -i /o nhl/2 „ 2 n+1/2 „ . 

_ rT.nd-1/2 , , *^2 . , / 3 R^ ^ 2 ^ sxn A0 

\,j.k - f\,3.k+ V ' -6^ 

. sin Ae, 

%>i,j,k - ~A0“> 


2 2 
+ (cos A0 - ^ ~ A6 sin_^. 

''gQ2^i,j,k 2 A0 ^ 


(47), 


Again, the order of accuracy is reduced slightly for 0^ larger 


than 45° 
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Due to the complex algegraic form of the generation terms, is 

just calculated from data at the intermediate time level. 

The above discussion briefly describes the numerical procedure. 
Additional detail can be found in the appendix at the end of this 
chapter. The numerical Integration is initiated by considering grid 
elements surrounding either pole. Since the area at the poles for 
iraput of mass is zero, the concentrations at the latitudes 10° from 
either pole do not directly depend on the concentrations at the poles. 

Representation of Sub-Grid Scale Motions ' 

Since the numerical approximation of the transport of mass attempts 
to represent a continuum by a discrete set of points, one must consider 
the effect of the small scale motions on the large scale convection. 

In any finite numerical model of turbulent flow processes, it is un- 
realistic to assume that motions of all scales can be simulated. 

Rather, one must be content to explicity evaluate the large scale 
convection and to parameterize the small scale mixing processes. 

This can be more readily appreciated by considering the appro- 
priately averaged species equation of continuity. Rather than time 
averaging the equations in the Reynolds sense as is customarily done, 
it is more useful to volume average over a given grid. Time averaging, 
in essence, represents all turbulent processes by the time averaged 
product of the deviation variables (e.g., u’v' , u’c’ , etc). However, 
in numerical treatment of the equations of change some of the motions 
that would be considered turbulent in the Reynolds sense may be 
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calculated explicity in time if the grid size is sufficiently small. 
Therefore, it is expected that the eddy diffusivity one uses to approx- 
imate the sub-grid scale motions depends on the grid size. 

More specifically, let the grid-scale averaging operator (repres- 
ented hy an overbar) be defined as follows: 

/ u((j),9,r,t)r^cos0 d^ d6 dr (48) 


c((f),6,r,t) r^cos0 d(f) d9 dr (49) 

r'^cos 0A^A9A r Ar A0 A(f) 

The grid— averaged variable (u or c) is now a continuous function of space 
and time. If one lets u’ and c’ represent deviations from local grid- 
volume means, it is easy to recognize that terms analogous to the con- 
ventionally time averaged ones are obtained. However, there is the 
Important difference between the two terms as to what portion of the 
turbulent processes do the two represent. 

There has been considerable attention, relative to solving the 
Navier-Stokes Equation, devoted to the concept of using a grid size 
dependent viscosity (sometimes referred to as an aritficial viscosity) 
to describe the micro-scale phenomena . Deardorff , 

Crowley Manabe , , Smagorinsky , Holloway, and Stone and others 
have employed such viscosities for fluid dynamics calculations. The 


u = i — / / 

r'^cos9A({)A0A r Ar A0 


or 


c = 


f f f 
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eddy viscosity can be based on either a local kinetic energy dissipation 

rate or an average kinetic energy dissipation rate. In the former case 

the resulting eddy viscosity is non-linear while in the latter case 

there is a single value for the entire space, 

(57) (59) 

Smagorinsky, Manabe, and Holloway \ Leith ^ , and Monin and 

Zilitinkevich^^^^ have used the well-known dimensional analysis 
approach of Kolmogorov for three-dimensional homogeneous turbulence 
to relate the eddy viscosity to the local rate of strain. Let be 
a dimensionless constant and let be the cube root of the product 
of the mesh spacing (for rectangular cartesian coordinates, = 

(A • A * A then 

xi X2 X3 


, 9u. 


9u. 

'1^ 


9u. 


1/2 


9x . 

X 


)] 


(50). 


In Equation (50) the Einstein convention for tensor notation is used for 

(51) 

simplicity, Deardorff's^ calculation for turbulent channel flow 
indicated that = 0,10 was optimum, 

Leith, however, argues that Equation (50) is of dubious validity 
in global numerical models since the horizontal grid is much larger 
than the vertical grid and thus is inconsistent with three dimensional 
isotropy. Instead, he treats the grid scale as being two dimensional 
and uses dimensional arguments based on the cascade of vorticity (co) 
in two dimensions to arrive at 
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e = a 2 ^^ IvwlA^ . (51) 

where = (A * A Crowley’ experiments on wind driven 

2 X2 

ocean circulations suggest that = 0,05. 

It is clear that Equatiors (50) and (51) lead to non-linear eddy vis- 
cosities which are local. Monin^^^^ has suggested that a linear viscosity 
can be estimated as 


£ ~ e 


1/3 

0 



(52), 


where is the kinetic energy dissipation averaged for the whole cal- 
culation space. For the entire atmosphere, ~ 5 erg/g sec. To intui- 
tively extend this concept one might consider e to be axis dependent and 
write 


£ 

X. 

1 




(53). 


While the simplicity of the linear viscosity is retained for computational 

purposes. Equation (53) does permit different scaling of the micro-scale 

(51) 

phenomena in the three directions. The results of Deardorff show 

4/3 

the constant of proportionality in Equation (53) to be (0.094) 

The discussion of the representation of sub-grid scale motions 
has thus far been restricted to momentum exchange. However, since our 
interests are in mass exchange the extension to eddy diffusivity must 
be made. The obvious appraoch is to assume that the turbulent Schmidt 


i 
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number (the ratio of the eddy viscosity to eddy diffusivity) is one. 

In addition, further detailing of approximations for the eddy diffusivity 

is hardly jusitifed at the present until more and better atmospheric 

turbulent exchange data are obtained. Based on these discussions, the 

following values for the eddy diffusivities would be consistent with 

grid spacings of 10“ x 10® x 2.5 km: £,=£„= lO^m^/sec and e = 10^ 

(p fct r 

m^/sec. These values are in agreement with the suggestions of Lilly 
and Deardorff and with the linear viscosity used by Crowley 
for ocean circulation calculations employing a smaller grid size, 

SUMMARY 

In this chapter, we have discussed the .development of a global 
transport/chemistry circulation model that accounts for natural and 
anthropogenic sources of ppllutants. Specifically the geochemically 
important methane - carbon monoxide system is described in which the 
oxidation of methane leads to the formation of carbon monoxide. The 
pseudo-steady state approximation is exploited for the reactive inter- 
mediate species so that the number of species continuity equations 
that must be solved is reduced to a minimum. By so doing three dim- 
ensionality of the solution can still be retained. 

The incorporation of the physico-chemical behavior of the species 
into the model is discussed in sufficient generality that the proce- 
dures described could be readily extended to other important systems. 

Such efforts should be encouraged. In the subsequent chapter, the 
results to date of computer simulations of the methane-carbon monoxide 


system are described 
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APPENDIX A 


DETAILS OF NUMERICAL SOLUTION 

In order to conveniently express Equation (40) we will define 


-oamea flux ter„s for ora These are 


,e 


Ft 
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(A-3). 


Let Xs y, and z represent the mole fraction of methane, formaldehyde, 
and carhon monoxide, respectively. The concentration of each species 
at the updated time can then be found by solving the following system 
of equations . 
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axe computed aud tlien Equation (A-4) is solved. Similar 

jK 

relationships can he written for the formaldehyde and carhon monoxide 
concentrations . 


Grid boxes that occur on a boundary of the region being considered 
can be treated using the integral flux method in a similar manner. As 
an example, consider the surface boxmdary where there is a grid box of 
dimension 2AcJ), 2A9 and Ar, By using the integral formulation, we 
can express the increment in the concentration of a particular species 
as 
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It should he noted that values of x, . should he regarded as 

i » 0 » 1 

averages over boxes centered on <j> 0^, ^3/2* j 2 

expressed in a central finite difference form, which employs stored 

values, as follows j 


jf 


i, j,2 


“ T~ [2 G. . - G. . (x” T )] (A-6) 

^ 3.»J»3 i,J»3 1,1,1 i,j,l 1,1,1^-^ 


Furthermore, the ground level flux terra must account for both a 
constant area source (Sg ) as well as an ocean or soil type source. 
Thus, 
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which to a fourth-order approximation can he written as 
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The convective fluxes can he conveniently expressed as 
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The vertical convective flux F: . „ can be obtained by Equation (44), 

1 j J 5(i 

and the generation term will be written as 
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Analogous to Equations CA~1) thru (A-eS) , define 
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With these equations, the concentration at the new time level for a 
boundary grid can be calculated as follows 
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Similar developments are made for the tropopause. 



43 


CHAPTER 2 

RESULTS OF PRELIMINARY SIMULATIONS 

The model described in the previous chapter was evaluated in static 
and preliminary dynamic transport simulations. Evaluations using the 
chemical-kinetic model in a static batch system analysis showed the chem- 
ical half lives of CH^ and CO, their non-chemical source strengths 
requisite with reasonable steady state concentrations, and the homo- 
geneous chemical source strength of CO by the gas phase oxidation of 
CH^ to be in the proper magnitude. In addition, analysis of the pseudo- 
steady state approximation for all of the""intermediates, including for- 
maldehyde, confirmed its utility and accuracy for simplifying the reaction 
scheme when one is not Interested in the very short term transient be- 
havior of the intermediates. 

Preliminary simulations using the transport/ chemistiry model were 
conducted, and the results of this early analysis are described. Chaining 
of the isopleths of CH^ and CO, due to the dominant westerly x^ind field, 
was observed. Longer term simulations are being planned. 

INTRODUCTION 

In the previous chapter the global transport/chemistry model for 
the CH^ - CO system was described. The physico-chemical considerations 
and mathematical development of the numerical model were presented in 
substantial detail. In this chapter, we will discuss the results of 
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computer simulations of the model and attempt to compare those results 
with some of the observations made on the species, especially carbon 
monoxide. 

Static and dynamic transport computer simulations were made. The 
chemical reaction model was first evaluated assuming a homogeneous 
reaction system of uniform temperature and pressure without considering 
any transport of the species. This was a preliminary analysis to judge 
the validity and evaluate certain aspects of the chemical reaction scheme. 

Of more significance is the combined transport/chemistry model simulation 
in which the troposphere’s variable properties and the Earth's variable 
surface are considered. The results of the static and dynamic model 
computer simulations will be discussed separately in the following sections. 

STATIC SIMULATIONS OF CHEMICAL REACTION MODEL 

The chemical kinetic model requires as input data the individual 
reaction rate constants and their temperature dependence, the third 
body concentrations, and the concentrations of the other chemical 

I 

entities that are considered time invariant. Rather than using these 
quantities as adjustable parameters, it was felt that these should be 
based on the best estimates from the literature. In Table I the 
reaction rate data and the appropriate references are listed. It should 
be noted that the temperature dependence is not known for all of the 
reactions. One should also note that the rate constants k, , k^, koj 
kg, and k^^^* required for the model when the pseudcv- steady state 
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TABLE I 

Reaction Rate Constants Employed in Model 


Reaction 

O^+hv ->■ 0(^D)+02 

0(^D)+H20 20H 

0(1d)-W^ -4- Of-M^ 

Cm)2+M^ Og+M^ 

CH^+OH CH 2 +H 2 O 

CH2O+ hv -»■ CO+H2 

CH 2 O+ hv + CHCH-H 

CH2O+OH + CHO+H2O 

CO+O+M. , CO„+M- , 
14 2 14 

CO+OH C02-+H 

Surface 

CO+N 2 O 


Rate Constant 

+1.05 X 10“5 e~°.48/cos0 ^^^-1 

3.0 X 10^^ m^/fcmole sec 

4.8 X 10^*^ m^/kmole sec 

3.0 X 10^ e m^/kmole^ sec 

T_10 -2500 /T 3., , 

2.0 X 10 e m /kmole sec 

+T cc -O.48/cos0 -1 

■1.65 X 10 e sec 

/o no”5 -O,48/cos0 -1 

'6.42 X 10 e sec 

, , mlO -460/T 3., - 

4.6 X 10 e m /kmole sec 

, , ,_11 -1750/T 6„ ,2 

3.6 X 10 e m' /kmole sec 

„ ^ -300/T 3,, , 

3.1 X 10 e m /kmole sec 

-6000/T -1 

24. e sec 


Reference 


2, 63, 64 

64, 65 
2, 64 

65, 66 

67, 68, 69 

70 

70 

64, 65, 67 
2, 65 
67 
2 


+Exponentlal term approximates the dependence of the reaction rate on the 
solar zenith angle; During the night, the exponential term was set to 
zero so that the photochemical reaction rates were zero. 
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approximation is used. The pairaraeteric data assumed for the other species 
are listed in Table II. 

The simplified reaction kinetic model was evaluated under static 
conditions to determine If it was generally consistent with the conclu- 
sions of other investigators. The reaction model was evaluated using 
CSMP on the IBM 370/165 system at the University of Kentucky. The 
Runge-Kutta fourth order variable step size integration routine was 
used. The temperature was specified to be 288°K and the water vapor 
concentration corresponded to approximately 35% relative humidity. 

The ozone mixing ratio was 0.02 ppm. 

The following five features of the reaction kinetic model were 
investigated; (a) analysis and validity of the pseudo-steady state 
approximation; (b) extension of the pseudo-steady state' approximation 
to include formaldehyde; (c) the chemical half lives of CH^ and CO; 

(d) the non-chemical source strengths of CH^ and CO requisite with 
reasonable steady state concentrations of CH^ and CO; and (e) the 
homogeneous chemical source strength of CO by gas phase oxidation of 
methane. Each of these will be discussed separately. 

Pseudo-Steady State Approximation 

It is virtually impossible due to computation time limiations 

to rigorously test the pseudo-steady state approximation by integrating 

the system of differential equations for total time periods on the order 
g 

of 19 seconds using a time step sufficiently small to ensure numerical 
stability and accuracy. However, an indirect validation can be achieved 

1 

!: 



47 


TABLE II 

Additional Parameter Values Employed in Model 


Parameter 


Temperature 


Pressure 


Solar Zenith Angle 


Assumed Value 
+ 288 ‘’K 

^1 atmosphere 

+ 43 ° 


0.02 ppm 

• —10 3 ■- 

.(8,4 X 10 kmoles/m at 

1 atmosphere and 288°K) 


H^O 




M, 


+2,5 X 10 ^ kmoles/m^ 

+Molar density of air 
0.042 kmoles/m^ 

+Molar density of air 

3 

0.042 kmoles/m 


^4 


+Molar density of air 

3 

0.042 kmoles/m 


Reference 


2, 64 


65 


65 


+These parameters were assumed constant in the static model, but were 
variable in the dynamic transport/chemistry model. 
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by evaluating the time required for the species to reach a relatively 

stable value and also to establish if this stable value corresponds to 

that determined when the pseudo-steady state approximation is employed. 

Several of the intermediate species were tested in this manner 

using the complete reaction mechanism without exploiting the pseudo-steady 

state approximation. Transport of the species was not important for 

this evaluation. The system of equations was integrated by evaluating 

the time step required to ensure numerical stability. In the initial 

-10 

stages, the time step was specified to be 10 seconds. When a species 

reached at least 95% of its pseudo-steady state value, the pseudo-steady 

state approximation was assumed for that species. In this manner, the 

/ 

time step could be increased slightly. All intermediate species were 
not tested since the computer time required would be impractical. 

However, sufficient data were obtained to rather conclusively show 
that the pseudo-steady state approximation is valid for this system. 

In Table III are listed the steady state values for the prescribed 
conditions and the time to reach 95% of the steady state value for 
three of the species. In the computer simulation, none of the species’ 
concentration exceeded the pseudo-steady state value. In addition, 
even those species that did not reach steady state in the allotted 
time appeared to be approaching the steady state value. Finally, it 
should be noted that the methane, formaldehyde, and carbon monoxide 
concentrations did not change during this time period. With this 
evidence, it is felt that the pseudo— steady state approximation is 
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Species 


0^(D) 


CH3O2 


TABLE III 

Test of Pseudo-Steady State Approximation 


Steady State 
Value 

2,44 X 10 ^^kmoles/m^ 
2.21 X 10 ^^kmoles/m^ 
5.81 X 10 ^^kmoles/m^ 
1.61 X 10 ^^kmoles/m^ 
6.70 X 10 ^*^kmoles/ni^ 
5.11 X 10 ^^kmoles/m^ 
7.10 X 10"^^kmoles/m^ 



4.05 X 10 ^ sec 
4.51 X 10 ^ sec 
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a valid modeling concept for this reaction scheme. 

Pseudo-Steady State Approximation for Foinxialdehyde 

With the success using the pseudo-steady state approximation for 

the intermediate species, this concept was also tested for formaldehyde 

since this species was not one of our principal interests. Long term 

3 

Integrations using integration time steps as large as 8.64 x 10 
seconds were made to compare the systems in which the pseudo-steady state 
approximation was either employed or not employed for the formaldehyde. 
Integrations for periods of 4.75 years showed differences in the con- 
centrations of less than 1%. In addition, the severity of numerical 
stability problems as evidenced by the required integration time step 
was much less when pseudo-steady state was assumed for formaldehyde. 

The reason for this is that even the three differential equations 
written for CH^, CH 2 O, and CO represent a stiff system, CH 2 O being the 
stiff variable. If the initial condition for CH 2 O does not closely 
correspond to the pseudo-steady state value consistent with the methane 
and carbon monoxide initial concentrations, small time steps must be 
used until the formaldehyde can adjust. On the other hand, this 
stiffness problem is circumvented if the pseudo-steady state approximation 
is employed for formaldehyde. Certainly one can not obtain the short 
term initial transient behavior (less than ten days or so) of formal- 
dehyde if the pseudo-steady state approximation is used, but the long 


term behavior is consistent with the transient model. 
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Chemical Half Lives, 

The chemical half lives must be at least. as long as the estimated 
overall half lives of CH^ and CO so that Inordinately high conversion 
rates are not simulated. There are various estimates for the overall 
residence times of methane and carbon monoxide. However, values of- 
2 *“ 3 years and 2-3 months, respectively, are not unlikely. For 
the purpose of testing this aspect of the chemical kinetic model, long 
term simulations up to 4000 days were conducted employing the pseudo- 
steady state approximation for all the species except methane and carbon 
monoxide. 

The results of those simulations are shown in Figure 3. The initial 

—8 3 

CH^ concentration was 6.3 x 10 kmoles/m Cl* 5 .ppm mixing ratio), and 

-9 3 

the initial CO concentration was 4.2 x 10 kmoles/m (0.1 ppm mixing 
ratio). The corresponding chemical half lives were 5.8 years and 1,1 
years for the CH^ and CO, respectively. These values are reasonable 
when one considers that sinks other than chemical reactions are also 
operating. 

Source Strengths of CH^ and CO 

If the troposphere can be considered a well mixed system, one 
can obtain estimates of the non-chemical reaction source strengths 
consistent with the chemical kinetic mechanism. This can lie accomplished 
by adding homogeneous sources of methane and carbon monoxide to Equa- 
tions (22) and (24), respectively, to account for all other sources of 
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these species. Then the magnitude of these source strengths is that 

which is consistent with steady state concentrations of CH^ and CO that 

one expects in the troposphere. The tropospherlcally averaged concen- 

_8 3 

trations of CH^ and CO are approximately 6.3 x 10 kmoles/m (1.5 ppm 

—9 3 

mixing ratio at 288®K)and 4.2 x 10 kmoles/m (0.1 ppm mixing ratio at 

288“K) . In Figure 4, the results of long term simulations where the 

—16 3 

homogeneous source strengths were fixed at 1.3 x 10 kmoles CH^/m 
-sec (320 X lO” tons CH^/year) and 1.5 x 10 kmoles CO/m sec 
(650 X 10^ tons CO/year) are shown. It is observed that with these 
values of the source strengths, the concentrations remain very nearly 
1.5 ppm and 0,1 ppm mixing ratio throughout the long term simulation. 
These source strengths would be changed slightly by assuming a different 
temperature and/or mixing ratio representative of the troposphere. 
However, these values are characteristic of the general magnitude that 
one anticipates for the source strengths of CH^ and the non-chemical 
reaction source strength of CO, 

CO Source Strength by Methane Oxidation 

An estimate for the source strength of carbon monoxide by the gas 
phase oxidation of methane can also b,e obtained. One procedure is by 
recognizing that, according to the simplified reaction scheme, all 
methane must ultimately be reacted to carbon monoxide. Thus, presuming 
typical methane and carbon monoxide concentrations and presuming that 
the* formaldehyde concentration can be specified by the pseudo-steady . 



Methane Concentration x 10® ( kmoles/m^ ) 
Carbon Monoxide Concentration x|0® (kmoles/m^) 



500 1000 1500 2000 2500 3000 3500 4000 

Time (Days) 

Figure 4 

Simulation for Source Strengths of Methane 'and Carbon Monoxide 
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state approximation, then the chemical source strength of CO determined 

■”16 3 

by Equation (22) is 1.3 x 10 kmoles/m"^ -sec. This corresponds to 

approximately 560 x 10^ tons CO/yr which is in reasonable agreement 

with other estimates. One can also note that at steady state, the 

homogeneous source of CH^ is totally converted to CO. As mentioned 
—16 3 

above, 1,3 x 10 kmoles CH^/m sec is also the best estimate and 
is in agreement with other work. 

Based on all of these evaluations, it can be stated that the 
simplified mechanism in combination with the pseudo-steady state 
approximation for all intermediates affords a plausible description of 
the CH^ - CO chemistry. With these significant simplifications the 
combined chemistry/transport model becomes much less time consuming 
on the computer. As a result of the pseudo-steady state approximation, 
it is only necessary to solve two species continuity equations, 

SIMULATIONS OF DYNAMIC TRANSPORT /CHEMISTRY MODEL 

The general approach in the simulation of the global transport/ 
chemistry model for the CH^-CO cycle consisted of the following steps; 
a. Initialize the CH^ and CO concentrations in the troposphere. 

In order to conserve computer time, the initial concentrations 
selected were approximately that expected in the atmosphere. 
This would not affect the final results but simply the com- 
puter time required to reach some regularly varying atmo- 
spheric state relative to the two species concentrations. 
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b. Distribute the sources and sinks of the various species on 
the Earth’s surface and at the tropopause consistent with the 
physico-chemical considerations. This involved proper inter-' 
pretatlon of oceans and lands as sources and/or sinks of the 
particular species. As a first approximation, the tropopause 
has been considered as a zero flux boundary. Thus, all 
pollutant generation and consumption is entirely within the 
troposphere . 

c. Solve the coupled unsteady state turbulent diffusion 
equatdons for CH, and CO with the boundary conditions 
established by b. Climatological data were used to estab- 
lish the wind field, temperature field, and water vapor 
field. . The coupling of the diffusion equations, of course, 
resulted from the gas phase reactions creating homogeneous 
generation terms. This, therefore, accounted for the chem- 
ical sources and sinks present, 

d. Continue the integration in time 

The inherent advantage to using this procedure is that one does 
not presuppose the atmospheric concentrations of the three pollutant 
species being studied. This, therefore, provides a meaningful test 
of the distribution of sources and sinks. This should not imply, 
however, that there are not uncertainties present. But these uncer- 
tainties are mostly associated with the strengths of the sources and 
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sinks, and a primary goal of the current research is to establish good 
estimates of these source and sink strengths. 

Model Parameters 

The basis for the model parameters was January climatological data 

f71 -j2') 

obtained from the National Weather Records Center ’ \ This 

included the horizontal wind field, temperature, and dew point at the 
surface, and the same data plus heights at pressure surfaces of 850, 700, 
500, 300, 200, and 100 mb. Dew points were available only through the 
500 mb surface. The computer model was not restricted to use of this 
climatology data. For example, the wind field could just as easily 
have been specified by a general circulation model. 

Since the global pollution transport model employed geometric 
height as the vertical independent variable rather than pressure, 
the data were converted to geometric altitude by linear interpola- 
tion. The vertical velocity was obtained from the horizontal wind 
field using the continuity equation. 

The chemical reaction rate constants were based on the tempera- 
ture field for those constants that showed a substantial temperature 
dependence and are listed in Table I in the preceding section. The 
chemical reaction model also required water vapor and ozone concen- 
trations as input. The water vapor concentration was established 
from the dew point data where possible. At levels of 5 km, 7.5 km, 

and -10 km, the water vapor concentration was assumed to be 2.4 x 10 ^ 

3 3 "6 3 

kmoles/m , 7.2 x 10 kmoles/m ,and 1.7 x 10 kmoles/m"^, in general 
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agreement with the U.S. Standard Atmosphere. The ozone mixing ratio 
was assumed to be constant at 0.02 ppm (v/v). 

The land and ocean areas were differentiated according to pollution 
source strength functions and their capacity to exchange the gaseous 
■pollutants with the air. Since the oceans apparently act as sources of 
both CH^ and GO, those locations assumed a supersaturation so as to 
describe the source. For methane, supersaturation corresponded to 1,8 
ppm for ocean areas and 1,0 ppm for land areas, whereas for carbon 
monoxide, these values were 3,5 ppm and 0.2 ppm, respectively. In those 
grids that had both ocean and land areas, a weighted average corres- 
ponding to the fraction of ocean and land area was used. Due to the 
limited knowledge of gas species interaction with the various soil 
types, there was no differentiation according to the type of land 
area. For example, desert areas were considered the same as forested 
regions in this initial simulation. Since the absorption coefficient, K, 
is dependent on the Henry’s Law constant which varies with temperature, 
the absorption coefficient for CH^ and CO were latitude dependent. 

The basis for the variation of H were zonally averaged sea surface 

(73) 

temperatures for January . The absorption coefficient was the 
same for land and ocean areas at the same -latitude. 

The anthropogenic sources were approximately .prorated according to 
the degree of urbanization. Therefore, the major amount of these 
sources were in the Northern Hemisphere. 
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It was necessary to consider the position of the sun to account 
for the photochemical reaction rate constants and the apparent diurnal 
nature of the supersaturation of CO in the oceans. The declination angle 
was set at -20“ consistent with about mid-January, and the initial 
time for the integration corresponded to 12:00 Greenwich mean time. 

Short Term Simulations 

Only one other simulation of this type appears to have been per- 

(74) 

formed. Kwok, Langlois and Ellefsen used a general circulation 

model and incorporated only anthropogenic, source estimates, with the 
atmosphere initially free of CO. In addition, their simulations did 
not include any possibility for homogeneous conversion of CH^ to GO. 

The present model has attempted to overcome these limitations. 

At the present^ only short term simulations have been performed, 
and longer term simulations are required. Typical results are shown 
in Figures 5 and 6. One of the important features, that of chaining 
of the Isopleths due to the dominant westerly wind field, can be 
observed being formed. Kwok, Langlois, and Ellefsen^^^^ also noted 
this with their CO transport model. 

Longer term simulations are being planned at this time. From 
those simulations, features such as the rate of interhemispheric trans- 
port, ground level source strengths, and the homogeneous conversion rate 
of CH^ to CO can be evaluated. These statistics will be used to validate 
and •further refine the model. 




^^:THA^£ 

time a 1.25 DAYS LE'vtL = 6250 METERS 


Figure 5 

Dynamic transport/chemistry model simulation of methane 
concentration. Units on isopleths are ppm x lo2. 



CARBON monoxide: 

TiM • 1.25 DAYS LEVEL » 6250 ME;TERS 


Figure ^ 

Dynamic transport/chemistry model simulation of carbon 
monoxide concentration. Units on isopleths are ppm x 103. 
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